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Abstract 

A method to compute the bound state eigenvalues and eigenfunctions of 
a Schrodinger equation or a spinless Salpeter equation with central interac- 
tion is presented. This method is the generalization to the three-dimensional 
case of the Fourier grid Hamiltonian method for one-dimensional Schrodinger 
equation. It requires only the evaluation of the potential at equally spaced 
grid points and yields the radial part of the eigenfunctions at the same grid 
points. It can be easily extended to the case of coupled channel equations and 
to the case of non-local interactions. 
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I. INTRODUCTION 



Numerous techniques have been developed to find the eigenvalues and eigenvectors of 
the Schrodinger and the spinless Salpeter equations. In particular, developments of the 
Hamiltonian in a convenient bases have been widely used (see for instance Refs. [llJ^). 
The accuracy of the solutions depends on two parameters: The size of the basis and a 
characteristic length which determines the range of the basis states. Upper bounds of the 
true eigenvalues are computed by diagonalizing the corresponding Hamiltonian matrix. The 
quality of the bounds increases with the size of the basis and, for a given number of basis 
states, there exist a characteristic length which minimize the value of a particular upper 
bound. 

In the case of a Schrodinger equation, other methods requiring only the evaluation of the 
potential at equally spaced grid points, yields directly the amplitude of the eigenfunctions at 
the same grid points I^Q. In particular, the Fourier grid Hamiltonian method appears 
very accurate and simple to handle. This method is variational and relies on the fact 
that the kinetic energy operator is best represented in momentum space, while the potential 
energy is generally given in coordinate space. 

In this paper, we show that this last method can be generalized to treat the semirelativis- 
tic kinetic energy operator, simply by developing the Fourier grid Hamiltonian method in 
the 3-dimensional space. Consequently, we propose to call our approach, the 3-dimensional 
Fourier grid Hamiltonian method. We focus our attention on the case of purely central local 
potential, but the method can also be applied if the potential is non-local, or if coupling 
exist between different channels. As explained below, the accuracy of the method depends 
on the number of grid points and on the maximal radial distance considered to integrate the 
eigenvalue equation. This last parameter is not easy to calculate without knowing a priori 
the wave function, so we propose an Ansatz to determine it. 

Our method is outlined in Sec. |I|, while Sec. |T| presents a convenient way to compute 
the domain on which the wave functions are calculated. Test applications of the method are 
described in Sec. and a brief summary is given in Sec. 0. 



II. METHOD 
A. Theory 

We assume that the Hamiltonian can be written as the sum of the kinetic energy T and 
a potential energy operator V . The eigenvalue equation for a stationary state is given by 



T + V 



m = E\m), (1) 



where T depends only on the square of the relative momentum p between the particles, V is 
a local interaction which depends on the relative distance, and E is the eigenenergy of the 
stationary state. This equation is a nonrelativistic Schrodinger equation if 

T = 1711+1712 + (2) 
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where mi and m2 are the masses of the particles and n is the reduced mass of the system 
(we use the natural units h = c = 1 throughout the text). Equation (|l]) is a spinless Salpeter 
equation if 



T = + ml + + ml (3) 
In configuration space, Eq. ([l|) is written 

J [{f\f\f') + {f\V\f')] (f'l^) df' = E (f 1^). (4) 

In the following, we only consider the case of a local central potential 

{f\V\r') = V{r)6{r—r') with r = |f |. (5) 

It is then useful to decompose the wave function into its central and orbital parts 

(r*|\E') = Ri{r) Yimir) with f = r/r. (6) 

To compute the non-local representation of the kinetic energy operator, we introduce the 
basis states {IkXu)}, which are eigenstates of the operator p. They are characterized by 
good orbital quantum numbers (A, u), obey the relation 

f{p^)\kXu) =T{e)\kXiy), (7) 

and satisfy the orthogonality relation 

{k'X'u'\kXu) = 6{k' - k) 6yx Su'u. (8) 

The representation of these states in the configuration space is given by 



{f\kXi^) = ^ — j^(kr)Yxuir), (9) 

where functions ji{kr) are spherical Bessel functions. Using the completeness relation of 
basis states {IkXu)} and Eq. (|^), we find 

{f\T\f') = / dk — Tik') Y: E Jxikr)j,ikr')Y,,if)Y:,if')- (10) 
Introducing the regularized function ui{r) = rRi{r), Eq. is written 

2 roc roo 

r dr'r'ui{r') dkk'^T{k'^) ji{kr) ji{kr') + V{r)ui{r) = Eui{r). (11) 



TT JO Jo 

This equation is the basis of the 3-dimensional Fourier grid Hamiltonian method. 



3 



B. Discretization 



We now replace the continuous variable r by a grid of discrete values defined by 

n = iA with i = 0, 1, . . . , N, (12) 

where A is the uniform spacing between the grid points. Regularity at origin imposes 
M/(ro = 0) = 0. For bound states, we have lim^^oo ^/('") = 0. Consequently, we choose to 
set ui{rN = NA) = 0. Actually, this last condition is not necessary but it does not spoil the 
accuracy of solutions. The normalization condition for the radial wave function is 

dr [uiir)f = 1. (13) 

The discretization of this integral on the grid gives 

A E i^iiri)? = 1- (14) 
1=1 

This corresponds to an integration by trapezoidal rule thanks to the choice of a vanishing 
radial wave function at r = tq and r = r^. 

The grid spacing A in the configuration space determines the grid spacing Ak in the 
momentum space. The maximum value of r considered being = NA, the wave function 
lives in a sphere of diameter 2rj\f in the configuration space. This length determines the 
longest wavelength Amax and therefore the smallest frequency Ak which appears in the k- 
space is 

We have now a grid in configuration space and a corresponding grid in momentum space 

ks = sAk = with s = 0, 1, . . . , A^. (16) 

If we note Vi = V{ri), the discretization procedure replaces the continuous Eq. (0) by an 
eigenvalue matrix problem 

N-l 

Y,H,,^] = en<P7 for 1 = 1,..., N-l, (17) 



where 



The (A^ — 1) eigenvalues e„ of Eq. (17) correspond approximately to the first (A^ — 1) eigen- 
values of Eq. (|TT|). In the case of a potential which possesses a continuum spectra, only 
eigenvalues below the dissociation energy are relevant. Other eigenvalues, which form a dis- 
crete spectrum of positive energies, are spurious and correspond to standing wave solutions 
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satisfying u{r) = at r = and r = A^A. The eigenvector 0" gives approximately the 
values of the radial part of the nth solution of Eq. ([TT|) at the grid points. The eigenvectors 
0" must be normalized according to Eq. ( p!^ in order that 0" ~ (r^). 

This method can also be used in the case of a non-local potential. If the interaction 
depends only on the radial variable, then the discretization of the action of the potential on 
the wave function gives 



7V-1 



dr' W{r, r') u{r') 



for 



(19) 



This corresponds also to an integration by trapezoidal rule thanks to the choice of a vanishing 
radial wave function at r = tq and r = r^- 

Coupled channels calculations can also be performed with this method. For instance, let 
us consider the coupled equations 

r ^(1) 10(1)) + 10(2)) = E|0(i)) 

\ |0(^)) + ^(2) 10(2)) = E|0(2)). 

The corresponding discretized equations are 



(20) 




(1) ^(1) 



f 
^(1) 



E, 
E. 



c(2) 



(21) 



(1 2) 

where Hlj ' and Wij are the 3-dimensional Fourier grid representation of the interaction 

are approximately the values of the radial part 
: iA for z = 1 - 1. 



operators i^(^'^) and W respectively, (pf'''^^ 
of the eigenstates |0(^'^)) at grid points rj = 



C. Relevance of discretization 

As shown in Sec. [11 A|, the 3-dimensional Fourier grid Hamiltonian method relies on the 



-xx' 
vr Jo 



following relation 

9 roo 

ji{kx) ji{kx') k'^dk = 6{x- x') . (22) 
The equivalent discrete orthogonality relation on our grid of points is 

s=l \ / \ / 

One can thus expect that A^^'''' = 6ij for all values of A^ and I. Actually, the situation is 
less favorable. As it is shown in the appendix, for Z = 0, we have 

A^f •'=°) = 5,,- VAT. (24) 



For / = 1, a|^'' 7^ 6ij, but we have verified numerically that 



'jj 7- 



hm A;f''=^) = 5,,. (25) 

For values of / larger than 1, formula ( pSj) is only approximately correct for small values 
of i and j. Consequently, the accuracy of this method becomes poorer when / increases; 
nevertheless for large enough number of grid points, very good results can be obtained. 
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III. DOMAIN OF INTEGRATION 



The accuracy of the eigenvalues and eigenfunctions depends on two parameters: The 
value of and the value of r^v- Obviously, for a given value of r^v, the accuracy increases 
with A^. A proper choice for the domain of integration is not evident. If r^r is too small, 
incorrect solutions will be found. If this parameter is too large, a great number of grid 
points will be necessary to obtain stable eigenvalues. In this section, we propose an Ansatz 
to compute a suitable value of r^r. The idea is to find the radial distance for which the 
radial part R{r) of the eigenfunction considered is such that 

r p/ ^1 < 26 

max [rK[r)\ 

where e is a number small enough to neglect the contribution of R{r) for values of r greater 
than r^. The eigenfunction considered being a priori unknown, we propose to use a trial 
wave function matching at best the true eigenfunction, at least for the large r behavior. The 
value of r satisfying the condition (^) for the trial wave function will be the value r^ used 
for the numerical computation. 

The first step is to find a potential Voo{r) which matches at best the potential V{r) for 
r ^ oo. In the following, we will consider three different types: 

Vooir) = K r^ with k> and p > 0, (27a) 
Voo(r) = with K>0 and < p < 1, (27b) 

V^{r) = -Vo e{a - r) with Vq > and a> 0. (27c) 

The second step is to choose a trial state |A) which depends on one parameter A, taken here 
as the inverse of a distance. This trial state and the eigenstate considered are characterized 
by similar behaviors for r ^ and r ^ oo. The best matching between this state and the 
trial state is obtained by means of the variational principle. The average value 

{X\H^\\) = {\\f + VUr)\\) (28) 

is then computed and the value of A is determined by the usual condition 

2<™=0. (29) 

In the case of the spinless Salpeter equation, the variational solution is computed using the 
fundamental inequality 



p + m?J <^{p)+m?. (30) 

The radial part R{r) of the trial state is then analyzed to find the value of r which satisfies 
the condition (pUj). 

We have remarked that with e = 10~^ it is possible to reach a relative accuracy better 
than 10~^ on eigenvalues, provided is large enough (A^ > 50 — 100). A relative accuracy 
on eigenvalues better than e can be achieved because the mean value of an observable is 
computed using the square of the function rR{r). 
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A. Ground states 



We first consider the case of ground states, tliat is to say states witliout vibrational 
excitation. In tlie case of a potential with a large r behavior given by Eq. ( |27a| ), we use 
harmonic oscillator wave functions as trial states. The radial part is given by 



R{r) 



r' exp 



(31) 



Using procedures ( pUj ) and Eq. ( PDD for the spinless Salpeter equation with potential 

( P7a|) , we find 



A 



p K- 



( 



\ 



+ 1) l^y(/ + |)A2 + mf ^{l + I) A 



2 + m2 



p+2 



(32) 



The corresponding relation for the case of a nonrelativistic kinematics is obtained with 
vanishing value for the parameter A in the right-hand side of the above formula. The 
reduced mass of the system appears naturally and the equation is no longer a transcendental 
equation. 

If the potential, at great distances, is similar to the potentials given by Eqs. ( |27bD -(27c), 
the trial states used are the bound state Coulomb wave functions. The radial part is written 



R{r) 



(2A) 



21+3 



■ exp(— Ar). 



\ r(2/ + 3) 

The variational calculation for the spinless Salpeter equation with potential ( |27b| ), gives 



(33) 



r(2/ + 3-p) 



+ 



r(2/ + 3) 

With the potential (|27c| ), we obtain 



2-p 



(34) 



A 



1 

2^ 



(2/ + 1) ln(2Aa) - In 



r(2/ + 3) 



A^ + mf V -^^ + "^2 



(35) 



Again, the corresponding relations for the case of a nonrelativistic kinematics is obtained 
with vanishing value for the parameter A under the square roots in the right-hand side of the 
above formulas. The reduced mass of the system appears naturally, but Eq. (^) remains a 
transcendental equation. 

Once A is found, it is easy to find r^. Let us introduce a dimensionless variable xn = ^r^. 
Using condition ( P^D with Eqs. (^Tj) and (^), Xa? is given by the transcendental equation 



(l + l) In 



X 



N 



l + l 



(36) 



with m = 2 in the case of Eq. (|3TD and m = 1 in the case of Eq. 
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B. Vibrational excited states 



When the eigenstate considered is characterized by a vibrational excitation v different 
from 0, we can use, in principle, the {v + l)th harmonic oscillator or Coulomb wave function 
as a trial wave function. But such a procedure makes analytical calculation of the optimal A 
much more complicated. One knows that the polynomial multiplying the exponential term 
in the {v + l)th wave function has degree {v + /) in the Coulomb case and {2v + /) in the 
harmonic oscillator case. So we can use a trial state with the value of / replaced by an 
effective orbital angular momentum /gfr which take into account the highest degree of the 
polynomial part of the radial trial state. We have verified that for potentials with large 
distance behavior of type ( |27a| ), it is a good approximation to take Ics = 2v + I. In the 
case of potentials with large distance behavior of types ( |27b| ) or ( |27cj ), it is better to use 

leS = V + I. 



IV. NUMERICAL IMPLEMENTATION 

We have tested the accuracy of our method with different models found in the literature 
HIJ^J^. In particular, we have find the same results as those of ref. in which a Schrodinger 
equation and a spinless Salpeter equation are used with a potential containing a Coulomb 
part and a linear part. In this section, we only present the results for a Schrodinger equation 
with a linear potential and for a spinless Salpeter equation with the Coulomb potential. 

In the model of Ref. [0], the masses of some meson states are simply given by a non- 
relativistic Hamiltonian with a confinement linear potential 

H = mi + 1712 + ^ + ar + C. (37) 

The regularized radial part M"(r) of the nth zero orbital angular momentum eigenfunction 
of this Hamiltonian can be written in terms of the Airy function 0] 



M'^(r) = (2/ia)i/6— ^=_=^, (38) 
//~Ai2(x) dx 



where Xn is the nth zero of the airy function. In Fig. [l|, we show the 8th S-wave eigenfunction 
of Hamiltonian (|37|) for parameters values: mi = 1712 = 0.300 GeV, a = 0.1677 GeV^ and 



C = —0.892 GeV, found in Ref. IQ (this corresponds to 7th excitation of the p-meson). On 
this figure, the exact function is obtained with formula ( pH]) and the numerical one has been 
computed with a value of calculated with the procedure described in Sec. |T| for e = 10~^ 
and with a number of grid points = 30. In these conditions, the 8th eigenvalue is found 
with a relative accuracy better than 10~^. This error can be reduced by a factor 10 or more 
by increasing A^. The numerical solution is indistinguishable from the analytical solution to 
the resolution of the figure. If the wave function must be used to compute mean values of 
observables, a greater number of points is obviously necessary. 

None analytical solution of the spinless Salpeter equation with Coulomb potential is 
known. But this equation has been extensively studied and it is possible to compare the 
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results of our method with results from other works. In Table |, we show some eigenvalues 
of the semirelativistic Hamiltonian 



with the parameter values: rrii = m2 = 1 GeV and k = 0.456. The accuracy of our results 
are similar of those of Refs. , even better for excited states found in Ref . (the purpose 
of the work in Ref. [@] was not to reach the greatest possible accuracy, but to demonstrate 
the feasibility of a method). We have remarked that a greater number of grid points is 
necessary for spinless Salpeter equation than for Schrodinger equation to reach a similar 
accuracy. 



The 3-dimensional Fourier grid Hamiltonian method, formulated and tested in this paper, 
appears as a convenient method to find the eigenvalues and the eigenvectors of a Schrodinger 
or a spinless Salpeter equation. It has the advantage of simplicity over all the other tech- 
niques. In particular, it requires only the evaluation of the potential at some grid points 
and not the calculation of matrix elements in a given basis. The method generates directly 
the values of the radial part of the wave function at grid points; they are not given as a 
linear combination of basis functions. Moreover, the extension of the method to the cases 
of non-local interaction or coupled channel equations is trivial. 

It is worth noting that the method based on the expansion of the wave function in basis 
functions can present some interesting features. In some cases, all the matrix elements can 
be generated from analytic expressions. Further, the size of the matrices required can be 
considerably smaller (about 20 x 20 or 40 x 40) 

The accuracy of the solutions of the 3-dimensional Fourier grid Hamiltonian method can 
easily be controlled since it depends only on two parameters: The number of grid points and 
the largest value of the radial distance considered to perform the calculation. A very good 
estimation of this last parameter can be easily determined by using the procedure described 
above, and the number of grid points can be automatically increased until a convergence is 
reached for the eigenvalues. The reliability of the method is also ensured by its variational 
character. 

The method involves the use of matrices of order {{N — 1) x (A^ ~ 1))? where N is the 
number of grid points. Generally, the most time consuming part of the method is the 
diagonalization of the Hamiltonian matrices. This is not a problem for modern computers, 
even for PC stations. Moreover, several powerful techniques for finding eigenvalues and/or 
eigenvectors exist and can be used at the best convenience. A demonstration program is 
available via anonymous FTP on umhsp02.umh.ac.be/pub/ftp_pnt/. 



We thank Prof. R. Ceuleneer, Dr F. Michel and Dr Y. Brihaye for useful discussions. 
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APPENDIX: ORTHOGONALITY CONDITION FOR S-WAVE STATE 

Using the development of spherical Bessel functions in terms of sine and cosine functions, 
we have 

{N,l=0) 2 . / TT .\ . / TT 

' — — > sm — St sm — . 



s=l 



Replacing the sine function in terms of exponential functions and using sin(O) = sin(7r) = 0, 
formula above becomes 



1 N-l 



s=0 

Distributing and using the well-known relation 



s=0 

one can obtain Eq. 



N-l 1 _ iirj 

E e'^'' = (A3) 
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TABLES 



TABLE L Energy eigenvalues of the spinless Salpeter equation with Coulomb potential 
V{r) = —k/t, for the parameter values mi = m2 = 1 GeV and k = 0.456. Our results, for three 

lo-^ 

IJ2l. 



values of with a value of rjv calculated with the procedure described in Sec. |I| for e 
are given with the upper bounds obtained by the variational methods described in Refs. 



State 



N = 100 



N = 200 



N = 300 



Ref. g\ 



Ref. § 



IS 
28 
38 
48 
IP 



1.9460 
1.9870 
1.9944 
1.9969 
1.9869 



1.9453 
1.9867 
1.9942 
1.9968 
1.9869 



1.9451 
1.9866 
1.9941 
1.9967 
1.9869 



1.9450 
1.9865 
1.9941 
1.9967 
1.9869 



1.9450 
1.9868 
2.0015 
2.0238 
1.9875 
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FIGURES 

FIG. 1. Comparison of exact (solid curve) and numerically computed (crosses surrounded by 
circles) eigenfunctions for the 7th excitation of the p-meson for the quark-antiquark Hamiltonian 
of Ref . 1^ . Our computation is carried out with = 30 and an integration domain determined by 
the procedure given in Sec. Ill for e = 10"*^. See the text for further details. 



13 




r (GeV^) 



